
options(max.print=5.5E8)
library(CDM)
##########################################################
# Cognitive Diagnostic Models in R                       #
##########################################################
pfad <- "E:\\DINA"
setwd(pfad)


#------------------------------------------------------------------
# Read in some example files
dinadat <- read.table( "data2017.csv" ,sep=",")
qmatrix <- read.table( "qmatrix.csv", sep="," )
content<-read.table("medical_items_2017.csv",header=T,sep=",")


# DINA estimation because rule = "DINA" is specified
d1 <- CDM::din( dinadat, q.matrix=qmatrix)
coef(d1)           

index<-d1$subj.pattern$pattern.index
cdm<-d1$pattern

plot(d1)
fitmod1 <-d1$item
fitmod2<-d2$item
write(t(fitmod1),file="E://DINA//fitdina.csv",ncolumns=6, sep= " ,")


################# Dimensionality Test  ###################################           
mod1 <- sirt::rasch.mml2(dinadat )
# estimate WLEs
wle1 <- sirt::wle.rasch(dinadat , b = mod1$item$b )$theta

# DETECT for content domains
detect<- sirt::conf.detect( data = dinadat , score = wle1  ,  itemcluster = content$content )
summary(detect)

mod2 <- sirt::rasch.mml2(dinadat)
wle2 <- sirt::wle.rasch(dinadat , b = mod2$item$b )$theta
# DETECT for Exploratory Perspective
           sirt::expl.detect( data=dinadat, score=wle2,
                nclusters=8, N.est=nrow(dinadat)  ) 
#############################################################
## Rasch Model ##############################################
#############################################################
# source routines
Rpath <- "E:\\DINA"
source( file.path( Rpath , "sirtr_0.9-06.R" ) )
setwd( Rpath )


library(ltm)
library(sm)
library(irtoys)
 library(sirt)


# Rasch model
ext1<-qmatrix$V1==1
ext2<-qmatrix$V2==1
ext3<-qmatrix$V3==1
ext4<-qmatrix$V4==1
ext5<-qmatrix$V5==1
ext6<-qmatrix$V6==1
ext7<-qmatrix$V7==1
ext8<-qmatrix$V8==1

c1<-subset(dinadat,select=c(which(c(ext1))))
c2<-subset(dinadat,select=c(which(c(ext2))))
c3<-subset(dinadat,select=c(which(c(ext3))))
c4<-subset(dinadat,select=c(which(c(ext4))))
c5<-subset(dinadat,select=c(which(c(ext5))))
c6<-subset(dinadat,select=c(which(c(ext6))))
c7<-subset(dinadat,select=c(which(c(ext7))))
c8<-subset(dinadat,select=c(which(c(ext8))))
                
##################################################################                
#########################  CTT model  ############################
##################################################################
ctt1<-apply(c1,1,mean)
ctt2<-apply(c2,1,mean)
ctt3<-apply(c3,1,mean)
ctt4<-apply(c4,1,mean)
ctt5<-apply(c5,1,mean)
ctt6<-apply(c6,1,mean)
ctt7<-apply(c7,1,mean)
ctt8<-apply(c8,1,mean)

par(mfrow=c(4,2))
hist(ctt1)
hist(ctt2)
hist(ctt3)
hist(ctt4)
hist(ctt5)
hist(ctt6)
hist(ctt7)
hist(ctt8)

CTT.data<-cbind(ctt1,ctt2,ctt3,ctt4,ctt5,ctt6,ctt7,ctt8)
CTT_mat<-matrix(NA,3259,8)
for(e in 1:3259) {
   for( f in 1:8) {
   if (CTT.data[e,f]>=0.5) CTT_mat[e,f]<-1 else CTT_mat[e,f]<-0}}



############  Rasch Person parameters #################################

 
map1 <- sirt::rasch.mml2( c1 , theta.k = seq(-3 , 3 , len=41) , conv1 = 0.001 , mitermax = 1)
fit1<-summary(map1)
r1<-map1$person$MAP
map2 <- sirt::rasch.mml2( c2 , theta.k = seq(-3 , 3 , len=41) , conv1 = 0.001 , mitermax = 1)
fit2<-summary(map2)
r2<-map2$person$MAP
map3 <- sirt::rasch.mml2( c3 , theta.k = seq(-3 , 3 , len=41) , conv1 = 0.001 , mitermax = 1)
fit3<-summary(map3)
r3<-map3$person$MAP
map4 <- sirt::rasch.mml2( c4 , theta.k = seq(-3 , 3 , len=41) , conv1 = 0.001 , mitermax = 1)
fit4<-summary(map4)
r4<-map4$person$MAP
map5 <- sirt::rasch.mml2( c5 , theta.k = seq(-3 , 3 , len=41) , conv1 = 0.001 , mitermax = 1)
fit5<-summary(map5)
r5<-map5$person$MAP
map6 <- sirt::rasch.mml2( c6 , theta.k = seq(-3 , 3 , len=41) , conv1 = 0.001 , mitermax = 1)
fit6<-summary(map6)
r6<-map6$person$MAP
map7 <- sirt::rasch.mml2( c7 , theta.k = seq(-3 , 3 , len=41) , conv1 = 0.001 , mitermax = 1)
fit7<-summary(map7)
r7<-map7$person$MAP
map8 <- sirt::rasch.mml2( c8 , theta.k = seq(-3 , 3 , len=41) , conv1 = 0.001 , mitermax = 1)
fit8<-summary(map8)
r8<-map8$person$MAP


par(mfrow=c(4,2))
hist(fit1[,3])
hist(fit2[,3])
hist(fit3[,3])
hist(fit4[,3])
hist(fit5[,3])
hist(fit6[,3])
hist(fit7[,3])
hist(fit8[,3])

write(t(fit1),file="E://DINA//fit1.csv",ncolumns=7, sep= " ,")
write(t(fit2),file="E://DINA//fit2.csv",ncolumns=7, sep= " ,")
write(t(fit3),file="E://DINA//fit3.csv",ncolumns=7, sep= " ,")
write(t(fit4),file="E://DINA//fit4.csv",ncolumns=7, sep= " ,")
write(t(fit5),file="E://DINA//fit5.csv",ncolumns=7, sep= ", ")
write(t(fit6),file="E://DINA//fit6.csv",ncolumns=7, sep= " ,")
write(t(fit7),file="E://DINA//fit7.csv",ncolumns=7, sep= ", ")
write(t(fit8),file="E://DINA//fit8.csv",ncolumns=7, sep= ", ")

########## Consistency among DINA, Rasch, CTT models  #############

cdm_dat<-data.frame(cdm[,2])

write(t(cdm_dat),file="E://DINA//cdm.txt",ncolumns=1, sep= " ")

examinees<-dinadat[,1]

  prop1<- matrix(1, nrow =length(examinees),ncol =1)
for (c in 1:length(examinees)) {
    prop=0
   for (d in 1:length(c1[1,])) {
     prop<-prop+(1/(1+exp((-1)*(r1[c]-map1$item$b[d]))))
     prop1[c]<-prop/length(c1)
     }}

   prop2<- matrix(1, nrow =length(examinees),ncol =1)
for (c in 1:length(examinees)) {
    prop=0
   for (d in 1:length(c2[1,])) {
     prop<-prop+(1/(1+exp((-1)*(r2[c]-map2$item$b[d]))))
     prop2[c]<-prop/45
     }}

    prop3<- matrix(1, nrow =length(examinees),ncol =1)
for (c in 1:length(examinees)) {
    prop=0
   for (d in 1:length(c3[1,])) {
     prop<-prop+(1/(1+exp((-1)*(r3[c]-map3$item$b[d]))))
     prop3[c]<-prop/45
     }}
     
    prop4<- matrix(1, nrow =length(examinees),ncol =1)
for (c in 1:length(examinees)) {
    prop=0
   for (d in 1:length(c4[1,])) {
     prop<-prop+(1/(1+exp((-1)*(r4[c]-map4$item$b[d]))))
     prop4[c]<-prop/25
     }}
     
  prop5<- matrix(1, nrow =length(examinees),ncol =1)   
for (c in 1:length(examinees)) {
    prop=0
   for (d in 1:length(c5[1,])) {
     prop<-prop+(1/(1+exp((-1)*(r5[c]-map5$item$b[d]))))
     prop5[c]<-prop/154
     }}     
     
 
  prop6<- matrix(1, nrow =length(examinees),ncol =1)   
for (c in 1:length(examinees)) {
    prop=0
   for (d in 1:length(c6[1,])) {
     prop<-prop+(1/(1+exp((-1)*(r6[c]-map6$item$b[d]))))
     prop6[c]<-prop/20
     }}     
     
  prop7<- matrix(1, nrow =length(examinees),ncol =1)   
for (c in 1:length(examinees)) {
    prop=0
   for (d in 1:length(c7[1,])) {
     prop<-prop+(1/(1+exp((-1)*(r7[c]-map7$item$b[d]))))
     prop7[c]<-prop/20
     }}     
     
  prop8<- matrix(1, nrow =length(examinees),ncol =1)   
for (c in 1:length(examinees)) {
    prop=0
   for (d in 1:length(c8[1,])) {
     prop<-prop+(1/(1+exp((-1)*(r8[c]-map8$item$b[d]))))
     prop8[c]<-prop/6
     }}     
      
     
rasch<-cbind(prop1,prop2,prop3,prop4,prop5,prop6,prop7,prop8)
rasch_dat<-matrix(NA,3259,8)

for(e in 1:3259) {
   for( f in 1:8) {
   if (rasch[e,f]>=0.5) rasch_dat[e,f]<-1 else rasch_dat[e,f]<-0}}
   

#write(t(rasch_dat),file="E://DINA//rasch.txt",ncolumns=1, sep= " ")

cdm_dat<-as.matrix(read.fwf(file="E://DINA//cdm.txt",width=c(1,1,1,1,1,1,1,1)))
#rasch_dat<-read.table("E://DINA//rasch.txt",header=FALSE)

dina_rasch<-matrix(NA,3259,8)   
for(e in 1:3259) {
   for( f in 1:8) {
   if (cdm_dat[e,f]==rasch_dat[e,f]) dina_rasch[e,f]<-1 else dina_rasch[e,f]<-0}}
 

CTT_rasch<-matrix(NA,3259,8)
for(e in 1:3259) {
   for( f in 1:8) {
   if (CTT_mat[e,f]==rasch_dat[e,f]) CTT_rasch[e,f]<-1 else CTT_rasch[e,f]<-0}}



CTT_dina<-matrix(NA,3259,8)
for(e in 1:3259) {
   for( f in 1:8) {
   if (CTT_mat[e,f]==cdm_dat[e,f]) CTT_dina[e,f]<-1 else CTT_dina[e,f]<-0}}

  
   
CC_dina_irt<-apply(dina_rasch,2,sum)/3259
CC_ctt_irt<-apply(CTT_rasch,2,sum)/3259
CC_dina_irt<-apply(CTT_dina,2,sum)/3259
   
   
cor(rasch)
cor(CTT.data)   
cor(cdm[,c(6:13)]) 
dcm<-cdm[,c(6:13)]

par(mfrow=c(4,2))
a1<-cor(r1,dcm[,1])
a2<-cor(r2,dcm[,2])
a3<-cor(r3,dcm[,3])
a4<-cor(r4,dcm[,4])
a5<-cor(r5,dcm[,5])
a6<-cor(r6,dcm[,6])
a7<-cor(r7,dcm[,7])
a8<-cor(r8,dcm[,8])  
DINA_IRT<-cbind(a1,a2,a3,a4,a5,a6,a7,a8)   

a11<-cor(CTT.data[,1],dcm[,1])
a22<-cor(CTT.data[,2],dcm[,2])
a33<-cor(CTT.data[,3],dcm[,3])
a44<-cor(CTT.data[,4],dcm[,4])
a55<-cor(CTT.data[,5],dcm[,5])
a66<-cor(CTT.data[,6],dcm[,6])
a77<-cor(CTT.data[,7],dcm[,7])
a88<-cor(CTT.data[,8],dcm[,8])  
DINA_CTT<-cbind(a11,a22,a33,a44,a55,a66,a77,a88)  

 
   
   
a111<-cor(CTT.data[,1],r1)
a222<-cor(CTT.data[,2],r2)
a333<-cor(CTT.data[,3],r3)
a444<-cor(CTT.data[,4],r4)
a555<-cor(CTT.data[,5],r5)
a666<-cor(CTT.data[,6],r6)
a777<-cor(CTT.data[,7],r7)
a888<-cor(CTT.data[,8],r8)  
IRT_CTT<-cbind(a111,a222,a333,a444,a555,a666,a777,a888)  

corr<-cbind(t(IRT_CTT),t(DINA_IRT),t(DINA_CTT))   
consistency<-cbind(CC_ctt_irt,CC_dina_irt, CC_dina_ctt)   
 